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A century on from its discovery, a complete fundamental understanding of su- 
perconductivity is still missing. Considerable research efforts are currently devoted 
to elucidating mechanisms by which pairs of electrons can bind together through 
the mediation of a boson field different than the one associated to the vibrations 
of a crystal lattice. PuCoGa 5 , a 5 /-electron heavy- fermion superconductor with a 
record critical temperature T c = 18.5 K, is one of the many compounds for which the 
short-range, isotropic attraction provided by simple electron-phonon coupling does 
not appear as an adequate glue for electron pairing. Magnetic, or virtual valence 
fluctuations may have an important role in the stabilization of the superconduct- 
ing ground state in PuCoGa 5 , but the specific nature of the coupling mechanism 
remains obscure. Here, we report the results of point-contact spectroscopy mea- 
surements in single crystals of PuCoGas. Andreev reflection structures are clearly 
observed in the low-temperature spectra, and unambiguously prove that the paired 
superconducting electrons have wavefunction with the rf-wave symmetry of a four- 
leaf clover. A straightforward analysis of the spectra provide the amplitude of the 
gap and its temperature dependence, A(T). We obtain A(T — > 0) = 5.1 ± 0.3 meV 
and a gap ratio, 2A/A;^T C = 6.5 ± 0.3, indicating that the compound is in the regime 
of strong electron-boson coupling. The gap value and its temperature dependence 
can be well reproduced within the Eliashberg theory for superconductivity if the 
spectral function of the mediating bosons has a spin-fluctuations-like shape, with a 
peak energy of 6.5 meV. Electronic structure calculations, combining the local den- 
sity approximation with an exact diagonalization of the Anderson impurity model, 
provide a hint about the possible origin of the fluctuations. 

PuCoGas, discovered in 2002, has the highest critical temperature (T c =18.5 K) among 
heavy-fermion superconductors Though the a-activity and the toxicity of Pu make it 
unsuitable for practical purposes, this material has being challenging the researchers since its 
discovery. As a matter of fact, it seems to elude any classification. The Sommerfeld coefficient 
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7, as measured by specific-heat experiments, is between 58 and 95 mJ/(mol K) |1|, [2|. The 
effective mass of electrons is about 1/5 of that measured for the isostructural unconventional 
superconductor CeCoIn 5 , suggesting a smaller degree of electronic correlation. 

The magnetic properties of PuCoGas and their possible relationship with the supercon- 
ducting pairing have been debated, mainly because of conflicting results about the magnitude 



2 



of the magnetic moment carried by the Pu ions. A Curie- Weiss susceptibility in the normal 
state was initially observed as expected for fluctuating local Pu 3+ moments. Successive 



/iSR studies showed no static, normal-state electronic magnetism in PuCoGa 5 [3j, |4j and polar- 
ized neutron diffraction experiments showed a small and temperature-independent microscopic 
magnetization dominated by the orbital moment [271 ] . 

The amplitude of the gap (or at least of the gap ratio 2A//cbT c ) has been estimated by 
different techniques. Using nuclear magnetic resonance (NMR) to measure the Knight shift in 
the superconducting state of a sample with T c =18.5 K, Curro et al. obtain 2A//cbT c = 8 (A ~ 
6.3 meV) 

From the theoretical and computational point of view, PuCoGas is a real challenge. Vari- 



8H13, 



23, 



291 ]. showing that the details of 



ous band structure calculations have been reported 
the Fermi surface (FS) and the value of the local magnetic moment depend on the theoreti- 
cal approach. In a local density /generalized gradient approximations (LDA/GGA), the FS is 



quasi-two dimensional 
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23], with at least two nearly-cylindrical sheets, one hole-like 



around the V = (0, 0, 0) point of the reciprocal space (k x , k y , k z ), and one electron-like around 
the M = (tt, 7r, 0) point. This bears a similarity to what more recently observed in Fe-based 
superconductors. The coupling mechanism that seems to account for the high-T c superconduc- 
tivity in iron pnictides (based on the nesting between hole and electron FS sheets through a 
vector associated with a peak in the spin susceptibility) had been proposed earlier for PuCoGa 5 



10]. The FS remains qualitatively unchanged with the addition of the Coulomb-?/ interaction 
29]. 

Several experimental facts 
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15Nl7l] suggest that in PuCoGas the electrons in 
the pairs have a mutual angular momentum i = 2, corresponding to a superconducting order 
parameter (OP) with d x 2_ y 2 symmetry, i. e. the gap in the single-particle excitation spectrum 
has a line of nodes intersecting the Fermi surface. A direct proof of the gap symmetry is, 
however, still missing. This is an important point, as the symmetry of the OP is closely 
related to the pairing mechanism. For instance, whilst isotropic electron-phonon attraction 
favors the formation of zero-angular-momentum pairs, with a spherically symmetric OP, d- 
wave symmetry is most easily realized if the pairing interaction is repulsive at short-range 
and anisotropic at larger distances, as the one provided by effective spin-spin couplings on the 
border of antiferromagnetism \\\ . 

A magnetic nature of the Cooper pairing mediator in PuCoGas has been discussed by several 
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authors (see Ref. 20J for a recent review). Flint et al. 21] considered virtual valence fluc- 
tuations of the magnetic Pu configurations, creating Kondo screening channels with different 
symmetry, and demonstrated that in a lattice of magnetic ions, exchanging spin with conduc- 
tion electrons in two different channels, a condensate of composite pairs between local moments 
and electrons is formed. These models must be reconciled with the temperature-independent 
susceptibility observed in the normal state 27]]. A phononic mechanism in the framework of 



the (i-wave Eliashberg theory has been discussed in Ref. [12J. Although able to reproduce a 
number of experimental observations, this model requires an electron-phonon (e-ph) coupling 
constant A that is much higher than the value experimentally deduced from the time relaxation 



of photo-induced quasiparticles 



22j (A = 0.2-0.26). 



We performed point-contact Andreev-reflection spectroscopy (PCARS) measurements in 
freshly annealed single crystals of 239 PuCoGas, grown by a flux method and characterized 
by x-ray diffraction, electrical resistivity, magnetization and specific heat measurements. The 
point contacts were made between a fresh, mirror-like surface of the crystal (just exposed by 
breaking the sample) and a thin Au wire (about 10 fim in diameter). The uneven broken sur- 
face on which the contact is made prevents a fine control of the direction of current injection 
with respect to the crystallographic axes. However, as shown below, this does not prevent the 
unambiguous determination of the amplitude and symmetry of the OP. 

Fig. [Upresents the temperature dependence of the raw conductance curves of a point contact 
whose normal-state resistance is 6.2 Q. All the curves but the lowest-temperature one are shifted 
downward for clarity. The top curve is measured at T = 1.8 K, that means about T c /10. This 
ensures that the gap extracted from it is well representative of the gap for T — > 0. The high- 
energy tails of the conductance curves are temperature independent, as shown in the inset, 
where the curves measured at T = 1.8 K and in the normal state are reported without vertical 
offset as an example. This demonstrates that the contact is in the perfectly ballistic regime, 
with no contribution from Maxwell terms in the contact resistance, i.e. there is no diffusion 
in the contact region and the maximum excess energy with which the electrons are injected 
into the superconductor is exactly eV. Meeting this condition is essential for energy-resolved 
spectroscopy to be possible. The absence of heating effects in the contact is also witnessed by 
the coincidence between the Andreev critical temperature T C A (i.e. the temperature at which the 
Andreev signal disappears and the normal-state conductance is recovered, here between 18.0 K 
and 18.5 K) and the bulk T c , as determined from resistivity measurements (T c = 18.1 K, see 
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Fig. 6 in the Supplementary Material). 

It is worth noticing that the point-contact curves in PuCoGas do not show the strong 
asymmetry observed for CeCoIn 5 below the temperature where a coherent heavy-fermion liquid 



develops 



35] . The asymmetric conductance in CeCoIn 5 has been explained in terms of a Fano 



resonance involving localized states near the interface and itinerant heavy electrons in the 
bulk ^M^. The absence of a strong, temperature dependent asymmetry in PuCoGas may thus 
confirm that superconductivity develops out of an incoherent metallic state j25| . 

In order to compare the experimental curves to the theoretical ones, a normalization is re- 
quired, i.e. a division by the conductance of the same point contact when the superconductor is 
in the normal state. Because of the extremely high value of the upper critical field in PuCoGas, 
this cannot be achieved by applying a magnetic field to suppress superconductivity at low tem- 
perature. However, owing to the moderate T c and the negligible temperature dependence of 
the normal state conductance (witnessed by the absence of any change in shape of the tails of 
the conductance curves at \eV\ > 15 meV) we can safely normalize the conductance curves at 
any T < T C A by the normal state conductance curve measured at (or just above) T C A . 

The result of this normalization is shown in Fig. [2] for two contacts made on different 
places of freshly-broken surfaces of the same sample. The two contacts have the same T C A 
but different normal-state resistance Rn- The curve shown in Fig. [2^ presents a very clear 
zero-bias conductance peak (ZBCP), whose amplitude (greater than 2) can only be explained 
by assuming a nodal OP with a change of sign at the Fermi surface. Indeed, the shape of 
the conductance curve is similar to those shown for a dimensionless barrier strength Z = 1 
in Fig. (Id) of the Supplementary Material, and clearly indicates that the OP has a 
c?-wave symmetry. The peak is ascribed to zero-energy Andreev bound states, arising from 
the constructive interference between electron-like and hole-like quasiparticles that feel order 

n 

parameters with different signs |3[. 

Another important point to notice is that, unlike in other heavy-fermion compounds like 
CeCoIns [36], the Andreev signal is very high in the superconducting phase; the curve is similar 
to the best ones ever observed in cuprates 291-132]. It thus seems that in PuCoGas there is no 
need to assume multichannel tunneling in the presence of localized states near the interface, as 
instead is the case in CeCoIn 5 j^. 

The conductance curve shown in Fig. [2b has a very different shape than the curve in Fig. 
reflecting a change from a = 7r/4toa = 7r/8of the angle between the direction of the 
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current injection (here in the k x -k y plane) and the k x axis. As a matter of fact, the curve 
in Fig. [2b looks similar to those shown in Fig. (lc) of the Supplementary Material, which 
have been obtained at T = for intermediate values of Z (0.5-0.7) and a misorientation angle 
a = 7r/8. 

Solid lines in Fig. [2] are the results of a quantitative spectral analysis using the 2-dimensional 
Blonder-Tinkham-Klapwijk (2D-BTK) formalism [l ^j]. The parameters of the model are A 
(here intended as the maximum amplitude of the gap), the barrier strength Z, the angle a, and 
a broadening parameter T that must usually be included in the model when fitting experimental 
curves. Here r < A so that it does not add any ambiguity to the determination of the gap 
amplitude. The best-fitting values of these parameters are indicated in the labels. The fit of 
the curves gives very similar values of the gap, and this is certainly the most important result, 
but also the barrier and broadening parameters are very similar. This means that the very 
different shape of the two curves almost completely depends on the angle a which is different 
in the two shown pictorially in the insets. 

Fig. shows an example of how the normalized conductance curves evolve on increasing 
temperature. The experimental data (symbols) are compared to the 2D-BTK fit (solid lines). 
The fit is rather good at any temperature; note that both Z and a are independent on T 
and were thus kept constant for all curves. We also kept T nearly constant so that the only 
parameter that varies significantly with temperature is the gap amplitude A. 

The temperature dependence of the gap as obtained from the fit of three different sets of 
conductance curves (in three different contacts) is shown in Fig. [3b. It is clear that the values 
are rather well reproducible; at low temperature, the gap values range between 4.85 and 5.30 
meV, corresponding to a gap ratio 2A//c B T c = 6.2 — 6.7. The vertical spread of gap values can 
be used to evaluate a posteriori the uncertainty on the gap itself; this is obviously much greater 
than the uncertainty arising from the fit of a single curve, which can be empirically determined 
as the range of gap values that allow an acceptable (i.e. within some confidence limit) fit of the 
conductance curve, when all the other parameters are varied as well. The spread of gap values 
is small at low temperature and maximum around 14 K. Although the general trend of the gap 
seems to be compatible with a BCS-like A(T) dependence (but of course, with a non-BCS gap 
ratio), this uncertainty does not allow discussing it in detail. 



The solid line in Fig. [3b is not a fit of the data, but is calcu 
coupling theory for superconductivity (known as Eliashberg theory 



ated within the strong- 
5]) by assuming that 
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the superconducting coupling is mediated by spin fluctuations. In particular, we used the 
typical spin-fluctuation spectrum with the form [8|: a 2 Fd{Vt) oc ^g^jm^n) -^here f2 is the 



energy peak, Q max is a cut off energy that we chose as being Q max = 4f2 [9|. The line in 
Fig. [3b is obtained by taking Q = 6.5 meV and using an electron-boson coupling constant 
A = 2.37. Details of the calculations are given as Supplementary Material. Assuming physically 
reasonable values of the Coulomb pseudopotential /i* (0 < /j* < 0.2) [5], the observed A(T) 
curve pose a stringent limit on the characteristic boson energy, namely 5.3 < Qq < 8 meV. 
This energy range is compatible with that determined by recent NMR measurements [13] for 
the characteristic spin-fluctuation energy in PuCoGas, that is 4 < < 8 meV. 

As a test of the reliability and generality of the results discussed so far, we also performed 
PCARS measurements in 242 PuCoGas crystals featuring a lower T c =14.5 K, due to the presence 
of Sb impurities (less than 1 %). The experimental results (shown as Supplementary Material) 
have been analyzed following the procedure described above, using the same spectral function 
and coupling constant as in Fig. [3b. We obtain a gap ratio 2 A/ ksT c = 6.2 ± 0.4, showing that 
A(T — >• 0) scales with T c . 

A scenario where magnetic fluctuations are responsible for the formation of the Cooper pairs 
in PuCoGas must be reconciled with the observed temperature-independent magnetic suscep- 
tibility [27] that points to vanishing local moments at the Pu sites. A plausible explanation is 
provided by electronic structure calculations combining the local density approximation (LDA) 
with the exact diagonalization (ED) [3^] of a discretized single-impurity Anderson model 16]. 
In this approach, the band structure obtained by the relativistic version of the full-potential lin- 
earized augmented plane wave method (FP-LAP W) 2l| is consistently extended to account for 
the full structure of the /-orbital atomic multiplets and their hybridization with the conduction 
bands [3]. Details on this procedure are given in the Supplementary Material. 

The resulting /-orbital density of states (DOS) is shown in Fig. 2^- Below the Fermi energy 
Ep, the DOS exhibits a three-peak structure that is typical for Pu and for a number of its 
compounds. Peaks characterized by the total moment j = 5/2 occur at Ep (peak labelled A) 
and at about 1 eV below Ep (peak C). The j = 7/2 peak B is located in between at approxi- 
mately 0.5 eV below Ep. The DOS is in a reasonably good agreement with the results of the 
non-crossing approximation reported in Ref. [261 ] . 

Fig-Hb shows the valence histogram calculated by projecting the ground state of the impurity 
model \Q) onto the Pu atomic eigenstates \m) that correspond to an integer /-shell occupation 
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n m . The plotted probabilities P m determine the /-orbital valence (n/) = 'Yl im P m n m = 5.3. 
The highest obtained probability is P5 = 0.63, followed by Pq = 0.33. In addition, there are 
small but non-zero probabilities P4 = 0.03 and P7 = 0.01. This result indicates a mixed- valence 
nature of the Pu / shell that consists of a mixture of a magnetic f 5 sextet and a non-magnetic 
/ 6 singlet, which is similar to the case of 5-Pu 43]. This mixed (also referred to as intermediate) 
valence is consistent with a large value of the 7-coefficient, as typically observed in the mixed- 
valence rare-earth-based materials. The same materials display a temperature-independent 
magnetic susceptibility at low temperatures if one of the involved atomic configurations is non- 



magnetic 



441 ]. which further strengthens the analogy between the mixed- valence rare-earth 



compounds and the /-electron physics in PuCoGa 5 . 

The Pu / shell carries a non-vanishing average moment as it fluctuates between the singlet 
and the sextet. The expectation values (X 2 ) = X(X + 1), where X = S, L and J, calculated 
for the / shell give S (5/) = 2.18, L (5/) = 4.05 and J (5/) = 2.43 for the spin, orbital and total 
moments. In the same time, the ground state \Q) of the entire impurity model is a singlet 
that corresponds to all angular momenta being equal zero (S — L — J — 0). Therefore, 
the fluctuations in the / shell are accompanied by canceling antiferromagnetic fluctuations 
in the conduction bands, which can be viewed as a manifestation of the Kondo physics. In 
addition, our band-structure calculations suggest an antiferromagnetic instability due to the 
presence of a Fermi-surface sheet with a negative second derivative of the Drude plasma energy 
(Supplementary Material). 

In conclusion, we have performed point-contact Andreev-reflection measurements in 
PuCoGas single crystals with different values of T c due to different degrees of disorder. We 
have shown that PCARS spectra unambiguously indicate that the OP has a c?-wave symmetry, 
consistent with indirect indications from NMR and /zSR measurements. The amplitude of the 
gap at T — > is A = 5.1 ±0.3, as determined from the fit of different curves, and corresponds to 
a gap ratio 2A/ksT c = 6.4 ± 0.4. In crystals with Sb impurities and a reduced T c = 14.5 K, the 
Andreev signal is greatly reduced, possibly because of a larger quasiparticle scattering in the 
superconducting bank arising from self-induced disorder and impurities. The gap scales with 
T c , giving the same gap ratio as in the highest-T c crystals, within the experimental uncertainty. 
In both cases, the temperature dependence of the gap is consistent with the predictions of the 
Eliashberg theory for strong-coupling superconductivity if spin fluctuations provide the medi- 
ating bosons. A characteristic boson energy Qq = 6.5 meV and a coupling constant A = 2.37 



S 



allow reproducing the low-temperature gap values. These results show that PuCoGas is an 
unconventional, rf-wave superconductor in the strong-coupling regime, and indicate spin fluctu- 
ations as the probable boson that mediates superconductivity. Such fluctuations would involve 
a time-dependent 5/ local moment dynamically compensated by a moment momentarily formed 
in the surrounding cloud of conduction electrons. 
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Supplementary material for "Strong-coupling d-wave superconductivity in 
PuCoGa 5 probed by point contact spectroscopy" 

I. POINT CONTACT ANDREEV REFLECTION SPECTROSCOPY 

Point contact spectroscopy is a simple but very powerful tool for the investigation of the 
superconducting order parameter (OP) that, in the past decades, has been successfully applied, 
among others, to cuprates, borocarbides, multiband superconductors and even to the recently 
discovered iron-based superconductors. The technique consists in creating a N-S contact be- 
tween a normal metal (N) and a superconductor (S), whose radius a is smaller than both the 
electronic mean-free path and the coherence length in S. In these conditions, an electron travels 
through the contact ballistically, i.e. without being diffused, so that if a voltage V is kept at 
the junction's end, it enters the superconductor with a maximum excess energy eV. When this 
energy is smaller than the gap in the superconductor, A, the electron cannot propagate in S 
as an electron-like quasiparticle (ELQ), because it finds no available states. Thus, it forms a 
Cooper pair and a hole is retro-reflected in N. 

If the gap is isotropic (Fig. [5^, top) and there is no potential barrier at the interface, this 
results in a doubling of the conductance for \ V\ < A/e (see Fig. [5b, top curve). If \eV\ > A, the 
conductance decreases again towards the value it would have if the superconductor was in the 
normal state. If a barrier is present, other phenomena take place that can give rise to a normal 
reflection of the incoming electron, and also to the transmission of hole-like quasiparticles (HLQ) 
in S. As a result, the conductance presents two maxima at approximately V = ±A/e and a 
zero-bias minimum (see Fig. [5b). 

If instead the energy gap has a dependence on the direction, i.e. the gap is anisotropic, 
electrons injected along different directions may experience different gap amplitudes. In this 
case, it is mandatory to integrate over all the possible directions along which the electrons 
can arrive to the interface. In particular, for a fi^.^-wave symmetry of the gap, A(#, <p) = 
Ao cos(2#) sin 2 (0), where 6 is the azimuthal angle in the k x ,k y plane of the reciprocal space 
and (f) the inclination angle (see Fig. [5^, bottom), the shape of the conductance curves does 
not depend only on the height of the potential barrier at the interface, but also on the angle 
between the direction of current injection (i.e. the normal to the interface) and the k x axis. 
This angle will be called a in the following. 
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The curves shown in Fig. [5]were calculated by using the so called 2D-BTK model [l-3] where 
the problem is reduced to a two-dimensional one and the Fermi surface (FS) is supposed to be 
perfectly cylindrical (with its axis parallel to the k z axis). This means that the dependence of 
the OP on cf) is disregarded. In this paper, based on the mostly 2D shape of the largest FS 
sheets and also for simplicity, we have always used this model. It can be shown |3J| that, with 
respect to a more refined 3D model, this approximation generally gives rise to an overestimation 
of the parameter Z, proportional to the height of the potential barrier at the interface, which 
is not relevant in our analysis. 

If a = 0, for any angle of incidence (8) of the incoming electron, ELQ and HLQ transmitted 
in S with angles +8 and —8 feel the same order parameter, in amplitude and sign. The 
conductance in this case is doubled only at zero bias where it shows a characteristic cusp. If 
there is no barrier, the same shape is obtained for any value of a: the conductance curve always 
looks like the top curve in Fig. [5b. If instead a barrier is present, for any a ^ some values 
of 8 exist for which HLQ and ELQ feel order parameters of opposite sign This gives rise 
to constructive interference between HLQ and ELQ that results in localized zero-energy states 
(Andreev bound states). These states manifest themselves in the conductance giving rise to a 
peak at zero bias (ZBCP). When a = tt/4 the current is injected along the nodal direction, all 
ELQ and HLQ interfere and the ZBCP is maximum. An example of calculated (normalized) 
conductance curves assuming a = tc/ 8 and a = tt/A is shown in Figs. [5b and[5bl for increasing 
values of Z. 



II. SOLUTION OF ELIASHBERG EQUATIONS 

In the imaginary- axis representation, the d-wave, one-band Eliashberg theory 0, 0] is for- 
mulated by the following equations for the gap A n (</>) = A(iu n ,<j)) and the renormalization 
functions Z n (4>) = Z (icu n , 0): 



n Z n (4>) = u n + nTJ2 r^\nm(<P, 00^(0') 



m 

TZ 



N n 

+ 7d c 2 + (A^) 2; (1) 



15 



ZMA n { ( f ) ) = nTj2 r% 

x[A nm (0, 00 -//*(0, <ff)&{u, c - K|)]iV£(0O, (2) 

where $ is the Heaviside function, w c is a cut-off energy for the Coulomb pseudopotential 
H*, and iV£(0) = A m (0)/v/^ + A^(0), A^(0) = w m / ^ + A m (0) 2 , and is an angular 
average of N% (00 over the Fermi surface. 

The parameter 7^ is proportional to the concentration of defects and c is a parameter related 
to the electron phase shift for scattering off an impurity [6]. The n th Matsubara frequency is 
defined as iu n = mT{2n — 1), and T is the temperature; A nm (0, 00 = \(ico m — icu n , 0, 00 is 
related to the electron-boson spectral function a 2 (Q)F(Q, 0, 00 through the following equation: 

\(iu m - iu) n ,(p,<p) = / n „ aii, (3) 

Jo {u m -u n y + iv 

where Q is the boson frequency. We made the usual lowest-order approximation that both 
the electron-boson spectral function and the Coulomb pseudopotential contain separate s- and 
d-wave contributions [7]: 

a 2 F(Vl, 0, 0') = a 2 F s (n) + 2a 2 F d (Q) cos (20) cos (20') (4) 

and 

/i*(0, 0') = fjt* + 2/i*(^) COS (20)c OS (20O (5) 

Then we assumed for simplicity that a 2 F ff (Q) = a 2 F d (Q), and that a 2 F d (Q) has the typical 

a]: a 2 F d (n) oc ""ogn^-n) where Qq ig the energy 

9| and w c = 3< >,„„.,. . 

We then looked for solutions of the equations having pure (i-wave symmetry for the gap 
function A(w,0O = A ( i(w)y / 2cos(20O and pure s-wave form for the renormalization function 
Z(u, 00 = Z s (u); the reason for this choice is that the only solution of the homogeneous integral 



shape of spin fluctuation spectral functions R ' • .»-/«'./(>» N "" !! " --'"-' :: - — 
peak, Q max is a cut off energy Q max = 4Q 



equation for Z d (u) is Z d (u) = 0, at least for reasonable values of A^ [l0|. Note that, if a d- 
wave symmetry for the gap function is assumed, the parameter //* does not enter into the two 
relevant Eliashberg equations. Therefore, although it is certainly larger than n d , it does not 
influence the solution. The free parameters are then: 7^, c, fi , fJ, a and A = J +o ° ma ^ d<yQ ^ dVL. 
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The parameters c and 7^ were estimated by reanalyzing the normalized local spin suscepti 



bility from the NMR experiment by Curro et al. 11]; in particular, it was necessary to assume 
7d = 0.25 meV and c = (unitary limit) to fit the shape of the curve 12]. The unitary limit 
is the only possible choice to reproduce the dependence of the critical temperature on disorder 



12J. As for /i d , it generally ranges between and 0.2 [5[. Finally, let us mention here that a 
range for VLq was recently provided by NMR [l3] measurements: 4 < VLq < 8 meV. 

The value of A was chosen so as to obtain the experimental critical temperature (T c = 18.5 
K) by solution of the imaginary-axis Eliashberg equations, and then the experimental low- 



temperature gap A by solution of the real-axis Eliashberg equations [14j for different values 
of fl and fi d . Fig. shows the values of the g function of VLq for different values of 

fi d . The corresponding values of A are shown in Fig. [6b. From the results of point-contact 
Andreev reflection measurement we know that, in PuCoGas single crystals with T c =18.5 K the 
low-temperature gap lies in the range 4.8-5.2 meV. As we can see from Fig. |6t, this corresponds 
to an interval of possible values of Qq that goes from 5.3 to about 8.9. From the intersection of 
this range with that given by NMR measurement [13J] one obtains 5.3 < Qq < 8 meV. Looking 
at Fig. [6b, it is clear that this corresponds to 2.2 < A < 3.7. 

To reproduce the temperature dependence of the gap, we chose for simplicity jj,* d — and 
adjusted the values of flo and A. We found that the experimental low-temperature gap and the 
T c can be obtained with Qq = 6.5 meV and A = 2.37. 

In order to obtain T c =14.5 K it is necessary to increase the scattering parameter to 7^ =1.6 
meV. In this way, we were able to reproduce the temperature dependence of the gap in the 
low-T c samples with the same values of S7 an d A mentioned above (see Fig. ITTb). This means 
that all the difference between the two sets of data can be ascribed to the greater amount of 
impurity scattering in the lower- T c crystals. 



III. ELECTRONIC STRUCTURE CALCULATIONS 



Following Ref. [15| we consider the multi-band Hubbard Hamiltonian H = H° + H mt , where 



H ° = E E H LnA^ = E E <, 72 (k) C t i (k)c 72 (k) 

k 71 ,72 



(6) 



hJ 71:72 

is the one-particle Hamiltonian found from ab initio electronic structure calculations of a peri- 
odic crystal. The indices i, j label lattice sites, 7 = (Ima) mark spinorbitals {</> 7 }, and k is a 
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vector from the first Brillouin zone. We assume that the electron-electron correlations between 
s, p, and d electrons are well approximated within DFT, and consider the correlations between 
the / electrons introducing the interaction Hamiltonian 

^ a,a' 

Rmt= 2^ ( m l ? m 2|V; ee |m3,m4)cL 1 a C l m2 a' C irn4 CT 'Cim3<x- (7) 

i mi,m2,TO3,m4 

The operator V ee represents an effective on-site Coulomb interaction [15| expressed in terms of 
the Slater integrals Fk and the spherical harmonics \lm). 

In what follows, the effects of the interaction Hamiltonian H mt on the electronic structure 
are modeled by a k- independent one-particle selfenergy E(z), which is constructed with the aid 
of an auxiliary impurity model describing the complete seven-orbital / shell including the full 
spherically symmetric Coulomb interaction, spin-orbit coupling (SOC), and crystal field (CF). 



The Hamiltonian of this multi-orbital impurity model can be written as [16 ] 

Himp — ^ ] [e ]mm'^ > krna^ > km'a' + ^ ] e ffmafma 

kmm' ma 
aa' 



mm' aa' 



+ J2{^mm'fLbkm'a'+h.C.) (8) 

kmm' 
aa' 

+ 2 ^ ] U mm i m ii m iiifl na f\ n , (y if m iii ( jif m ii a ^ 



mm'm" 
m"'aa' 



where f^ creates an electron in the / shell and creates an electron in the "bath" which 
consists of those host-band states that hybridize with the impurity / shell. The impurity- 
level position Ef and the bath energies e k are measured from the chemical potential /i. The 
parameters £ and Acf specify the strength of the SOC and the size of the CF at the impurity. 
The parameter matrices V k describe the hybridization between the / states and the bath 
orbitals at the energies e k . 



The band Lanczos method 
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18] is employed to find the lowest-lying eigenstates of the 
many-body Hamiltonian H imp and to calculate the one-particle Green's function in the subspace 
of the / orbitals [Gi mp {z)}^ n ^, at low temperature {ksT = 1/500 eV). The sought for selfenergy 
[£(2)]^^/ is then straightforwardly obtained from the inverse of the Green's function matrix 
[G\mv{ z )\mm' ■ Note that for the Lanczos method to be applicable, the continuum of the bath 
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states is discretized. In particular, the calculations presented here utilize only a single value 
for the index k. 

Once the selfenergy is known, the local Green's function G(z) for the electrons in the solid 
is calculated as 

[G(z)}^\ 2 = [G LDA (^)]- 1 72 - Ae5 7l72 - [E(z)] 7172 , (9) 

where Ae accounts for the difference between the impurity and the lattice chemical potentials, 
and Glba(z) is the LDA Green's function, 



[Glda(2)] 7i72 = tt- I d 3 k [z + n - F LDA (k)] 1 

VBZ JBZ 7 1 



(10) 



For the charge-density self-consistency, we employ the so-called local density matrix approxi- 



mation (LDMA) 19]. In the LDMA, the occupation matrix 



i r F 

-lm dz [G{z)\ lll2 (11) 

^ J-oo 



is self-consistently evaluated with the aid of the local Green's function G(z) from Eq. (Q. The 
matrix n 1112 is then used to construct an effective LDA+U potential Vu which is inserted into 
Kohn-Sham-like equations, 

[-V 2 + Vlda(i") +V n + £(1 • s)]<(r) = «(r). (12) 

These equations are iteratively solved until self-consistency over the charge density is reached. 
In each iteration, a new Green's function Gxda(z) and a new value of the /-shell occupation are 
obtained from the solution of Eq. (|T2l) . Subsequently, a new selfenergy £(z) corresponding to 
the updated /-shell occupation is constructed. Finally, the next iteration is started by inserting 
the new Glda(^) and into Eq. 0. After the iterations are completed, the imaginary part 
of the local Green's function G(z) provides a means to estimate the valence-band photoemission 
(PE) spectra. 



In the calculations we used the in-house implementation 2Cj, 1211] of the full-potential lin- 



earized augmented plane wave (FP-LAPW) method 22J. This FP-LAPW version includes all 
relativistic effects: scalar-relativistic and spin-orbit coupling. The calculations were carried out 



assuming a paramagnetic state with the LDA optimized 



231 ] crystal structure parameters of 



PuCoGa 5 . The Slater integrals were chosen as F = 4.0 eV, and F 2 = 7.76 eV, F 4 = 5.05 eV, 



and Fq = 3.07 eV 24]. They corresponds to commonly accepted values for Coulomb U = 4.0 
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eV and exchange J = 0.64 eV. The SOC parameter was determined from LDA calculations as 
£ = 0.28 meV. The crystal-field effects were found to be negligible and we set Acf = 0. 

The first and fourth terms in the impurity model, Eq. (jSJ), are to a good approximation 
diagonal in the {j, j z } representation, so that we only need to specify one bath state (six orbitals) 
with e^Z^ 2 and V k ~^ 2 , and another bath state (eight orbitals) with e^zJ/ 2 and V^~i ■ In 
order to determine these bath parameters, we assume that LDA represents the non-interacting 
model for PuCoGas, that is, we associate the LDA solution with the Hamiltonian of Eq. (jSJ) 
where we set the Coulomb interaction matrix U m m'm"m"' = 0. The numerical values of the 
bath parameters are then found as follows: We calculate the LDA hybridization function A = 
ImTr[GLDA( e + ^)] f° r 3 = ^/2 (Nf = 6) and for j = 7/2 (Nf = 8), which is shown in 
Fig. [71 Assuming that the most important hybridization is the one occurring in the vicinity 

C fey 7/0 ■ 2 J 

of Ep, we obtain V k L{ from the relation J2 k \ V£\ 5(e 3 k — e) = A J (e) [25] averaged over the 
energy interval e G (Ep — 0.5 eV, Ep + 0.5 eV). The bath-state energies e^ 2 ^ 2 are adjusted to 
approximately reproduce the LDA 5 /-states occupations n! 2 and n 7 / 2 . 



TABLE I: 5 /-states occupations rif, 


5/2 j 7/2 

nJ and nJ , 


el /2 ' 7/2 (eV), 


^5/2,7/2 


(eV) for PuCoGa 5 . 


5/2 7/2 

n f nJ nj 


5/2 
e l 


^5/2 


7/2 


vl 12 


5.14 4.38 0.76 


0.25 


0.29 


-0.07 


0.34 



In what follows, we have adopted the fully-localized (or atomic-like) limit (FLL) for the 
double-counting term V^ c = U(jif — l/2) — J(nf — l)/2 entering the definition of the LDA+f/ po- 
tential Vjj. In the FP-LAPW calculations we set the radii of the atomic spheres to 3.1 a.u. (Pu), 
2.3 a.u. (Co) and 2.3 a.u. (Ga). The parameter R x K max = 10.2 determined the basis set size, 
and the Brillouin zone (BZ) sampling was performed with 1152 k points. The self-consistent 
procedure defined by Eqs. fl9|)- f|T2|) is repeated until the convergence of the /-manifold occupa- 
tion n/ is better than 0.01. 

The resulting /-orbital DOS (or, the one-particle spectral function) is shown in Fig. 4a of 



the main text. The DOS below Ep is in reasonably good agreement with Ref. 26[. Above Ep, 
the DOS obtained here differs somewhat from the DOS reported in Ref. {26 1 for a reason which 
is unclear at the moment. We estimate the electronic specific heat coefficient using the formula 
7 = ^k 2 B Tr[N(Ep)(l — d ^^ )U=o] an d we get ~ 40 mJ K -2 mol -1 which is smaller than the 
experimentally estimated value of 80-100 mJ K -2 mol -1 . Note that a possible enhancement of 

20 



7 due to the electron-phonon interaction is not taken into account. 

Now we discuss the magnetic properties. The ground state of the cluster formed by the 
/ shell and the bath, Eq. (jSJ), is a singlet (S = L = J = 0) that includes (rif) = 5.30 in the 
/ shell, and (nb at h) = 8.70 in the bath states. The spin S = 2.18, orbital L = 4.05 and total 
J = 2.43 moments are calculated for the / shell from the expectation values (X 2 ) = X(X + 1), 
X = S,L, J. The individual components of the moments vanish, (S z ) = (L z ) = 0, unless the 
symmetry is broken by an external magnetic field. We carry out calculations with a small 
external magnetic field added to Eq. ([8]) in order to estimate the local /-shell susceptibility. 
These calculations yield ~ 0.3 x 10~ 3 /ie/T which is somewhat smaller than the experimental 
value 0.74 x 10" 3 fi B /T |27|. 

Qualitatively, the non-magnetic character of the composite electron cluster is linked to the 
mixed-valence nature of the Pu atom. The fluctuating local moment is compensated by corre- 
sponding dynamical fluctuations of the bath moment. The magnetic susceptibility x ~ (T+T f ) 
is anticipated, similar to a Kondo-singlet state, which remains constant for T « Tf c , as ob- 
served experimentally [27]. Note that neither the admixture of atomic-like / 5 and f 6 multiplets 
nor the Kondo-like screening of the f 5 state localized moment 28J by the bath alone are suffi- 
cient to form a singlet ground state of the composite cluster formed by the f-shell and the bath. 
When the hybridization is set to zero in Eq. ([8]) (LDA+Hubbard-I approximation 29] ) the 
resulting mixed- valence solution with (rif) = 5.3 is not a singlet and the f-shell stays magnetic. 
On the other hand, when we force f 5 solution in Eq. ([8]) and keep hybridization unchanged, 
the ground state of the composite electron cluster is not in a singlet either. In this sense, the 
admixture of the non-magnetic f 6 to the magnetic f 5 state enhances the dynamical screening 
of the local moment by the bath. 

Comparison of theoretical and experimental PE spectra represents an important criterion 
of soundness of the electronic structure calculations. The use of the single-particle LDA DOS 
to compare with the PE spectrum is justified only in case of weak electron correlations among 
the 5/ states. For the 5/ states at the borderline between the localized and itinerant behavior 
neither the LDA nor the static-mean-field LDA+U theories are sufficient to accurately describe 
the PE since these approximations entirely miss the experimentally observed atomic-like multi- 
plets. In Fig. [8]we show the total DOS for PuCoGas as provided by our LDA+ED calculations. 



Overall, the agreement with the experimental PE data 
calculations reported earlier 29]. 



30j is improved over LDA+Hubbard-I 
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Now we discuss the static band-structure resulting from the solutions of Eq. (fl2~|) . In fact, 
as it can be seen from Fig. [9] (upper panel), near Ep it is very similar to the static band 



structure obtained from LDA+HIA calculations 
surfaces (FS) from Ref. 



29| . The corresponding quasi-particle Fermi 



are shown in the lower panel of Fig. |HJ There are four sheets (I-IV) 
composing the FS: sheets I and II are fairly three-dimensional and sheets III and IV are two- 
dimensional. Close similarities are revealed between these sheets and the Fermi surfaces from 
the previous LDA 23] and AMF-LDA+U calculations. As it was already noticed in Ref. 29], 



FS calculated with LDMA and LDA have very similar geometry. 



Now we calculate the Drude plasma energy, Vl p = J + + where the individual 
directional components Q{ are given by 

^ = -f~J d'kv^EM-Ep), (13) 
Z7r Jbz 

which leads to numerical value of Q p = 3.8 eV, as in previous LDA+HIA calculations 29]. 
Eventually, it leads to the same estimate for the electron-phonon coupling A tr ^2.5 from the 
Bloch-Griineisen transport theory. 

In addition, we can estimate the type of the magnetic instability, ferro- or antiferromagnetic, 
for the band structure in Fig. [9) The paramagnetic susceptibility XoiQ^) m the limit of small 
q and u can be expanded [31] as Xo(Qi u ) = N(Ep) — aq 2 + lb-. In a simplified form [32], the 
coefficient a is proportional to a second derivative d \e=e f of the Drude plasma energy. The 
sign of this derivative indicates the type of magnetic instability. A positive sign corresponds to 
q = ferromagnetic fluctuations, whereas a negative sign to q ^ antiferromagnetic (or spin- 
spiral) fluctuations. In Table II we show the second derivative of Qp(E) for all four FS bands 
calculated from cubic polynomial fits. Our estimate points to FS-II as a suitable candidate for 
q 7^ antiferromagnetic instability. The electrons at FS-II, at least in principle, can be paired 
by anti-ferromagnetic fluctuations. 

It is also to note that the shape of FS shown in the lower panel of Fig. [9] suggests a possibility 
of nesting between the hole-like FS-1 around the T = (0, 0, 0) point and electron-like FS-3 
around the M = (tt, it, 0) point. This large q ~ (tt, tt, 0) nesting-like feature can lead to a 
peak in the spin susceptibility, and promote the anti-ferromagnetic fluctuations to develop. 
This suggest similarities between the coupling mechanisms for the high-T c superconductivity 
in PuCoGas and Fe-based superconductors [33 1. 
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TABLE II: The total BOS{E F ) in eV"\ Drude plasma energy (eV), and its second derivative over 
energy at Ep 

FS-I FS-II FS-III FS-IV 
DOS(£ F ) (1/eV) 0.7 1.3 0.8 0.5 
Q P (E F ) (eV) 1.1 2.7 1.9 1.5 
^"Iff 49.71 -221.14 72.0 193.7 

IV. CHARACTERIZATION OF THE INVESTIGATED SAMPLES 

The samples used for these experiments have been characterized by magnetic susceptibility 
and electrical resistivity measurements, using commercial Quantum Design platforms (MPMS 
7-T SQUID, and PPMS-9T). 239 PuCoGas single crystals have been submitted to a thermal 
treatment to anneal the self-radiation damage. Fig. [TO] shows the temperature dependence 
of the electrical resistivity measured for one of these crystals immediately after a thermal 
treatment. All the samples exhibited a critical temperature very close to the optimal value T c 
= 18.5 K. Point Contact Andreev Reflection Spectroscopy (PCARS) measurements were started 
within one day from the thermal treatment. PCARS measurements have also been performed 
on 242 PuCoGa 5 crystals featuring a lower T c =14.5 K, due to the presence of Sb impurities (less 
than 1 %). The results obtained for these samples are presented in the following section. 



V. PCARS MEASUREMENTS ON SAMPLES WITH T c =14.5 K 

The point contacts were made on the side surface of the crystals so as to inject the current 
mainly along the ah planes. Unlike in the purest samples, here we used a small spot of Ag 
conducting paste between the Au wire and the sample to act as the N electrode and also to 
mechanically stabilize the contacts. Fig. fTTh shows a series of conductance curves measured 
as a function of temperature in a contact with Rn = Q.9Q. The curves are vertically offset for 
clarity; as in the cases discussed for the sample with T c = 18.5 K, there is no shift of the tails 
on increasing temperature. 

Fig. [TTb reports the normalized low-temperature conductance of the same contact (top) 
and that of a different contact with a higher Rn =18.8 Q (bottom). The experimental data 
(symbols) are compared against the relevant 2D-BTK fit (lines). It is worth noting that these 
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curves (and all the curves obtained in crystals with reduced T^j show a small Andreev signal 



whose amplitude is similar to that observed in CeCoIns 



34, 



351 ] . It seems however logical, in 



this case, to ascribe the small signal to the greater amount of disorder in the samples (and thus 



also on the sur 
as proposed in 



ace) with respect to the purest crystals, rather than to intrinsic phenomena 



361 ] . The small amplitude and the consequent absence of very clear structures 



result in a greater uncertainty on the fitting parameters, since different sets of parameters can 
give equally good fits. Moreover, the fit requires a large T (although still smaller than A). As 
an alternative, and according to 36|, one could assume that a fraction of the injected current 
tunnel into a non-superconducting band or set of states, and does not contribute to Andreev 
reflection. In the latter case, the values of T would be smaller but an additional parameter 
would be necessary. 

Fig. [Tib shows an example of temperature dependence of the gap extracted from the fit 
(symbols). The line is the curve calculated within the Eliashberg theory, by using the same 
spectral function and coupling constant as for the sample with T c = 18.5 K. The different T c 
and the different gap amplitude arise only from the larger value of the quasiparticle scattering 
rate 7^ included in the Eliashberg equations to account for the disorder. Incidentally, while 
T c = 18.5 K is obtained by using 7^ = 0.25 meV, the reduced T c = 14.5 K of this particular 
contact requires 7^ = 1.6 meV. Note that also the gap ratio 2A/A;#T C = 6.2 ± 0.4 is compatible 
with that obtained in the purest crystals. 
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FIG. 1: Temperature dependence of the raw conductance curve of a Au/PuCoGas point contact with 
Rn = 6.2 All the curves apart from the top one are vertically shifted for clarity. The inset shows 
the curve at 1.8 K and the curve at 20.02 K (and thus in the normal state) without any shift. The 
superposition of the tails is one of the indicators that the contact is in the ballistic regime. 
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FIG. 2: Two examples of low-temperature conductance curves in PuCoGas after normalization (sym- 
bols) and the relevant 2D-BTK fit assuming a d-wave gap (lines). The values of the parameters are 
indicated in labels, as well as the normal-state resistance Rn of the corresponding Au/PuCoGas point 
contact. A is the gap at T = 0, T is a spectral broadening parameter, Z measures the barrier strength 
at the interface, and a is the angle between the direction of current injection and the k x axis. The 
different shape of the two curves simply arises from the different direction of the current injected with 
respect to the lobes of the ci-wave symmetry (insets). 
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FIG. 3: (a) Temperature dependence of the normalized conductance curves of a contact with Rn = 6.0 
(symbols) and their fit with the 2D-BTK model (lines). The curves are vertically offset for clarity, 
(b) Temperature dependence of the gap extracted from the fit of the conductance curves of three 
different contacts (symbols). The solid line is the A(T) function calculated within the Eliashberg 
theory by assuming a spin-fluctuation spectrum peaked at f^o = 6.5 meV and a coupling constant 
A = 2.37. 
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FIG. 4: (a) /-electron density of states (DOS) for the Pu atom in PuCoGas: j = 5/2, 7/2 projected, 
and total; (b) Valence histogram obtained by the projection of the LDA+ED solution onto Pu atomic 
eigenfunctions. 
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FIG. 5: (a) Amplitude of the order parameter in spherical coordinates in the reciprocal space. The 
cases of a s-wave (top) and of a cf-wave (bottom) order parameter are shown, (b): Theoretical 
normalized conductance curves calculated at T = for a s-wave gap, and different values of the 
dimensionless barrier strength Z. (c,d) Normalized conductance curves at T=0 calculated with the 
2D-BTK model for a d-wave gap and assuming a = n/8 and a = 7r/4, respectively; a is the angle 
between the direction of the injected current and the k x axis. 
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FIG. 6: (a) The effect of Slo ° n A for different values of the Coulomb pseudopotential fj,*,- The 
dashed region indicates the range of gap values extracted from the point contact Andreev reflection 
measurements. The inset shows the shape of the electron-spin-fluctuation spectrum, (b) The effect of 

Oo on the coupling constant A for different values of the Coulomb pseudopotential. 
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FIG. 8: (Color online) Total DOS of PuCoGa 5 provided by LDA+ED with U = 4.0 eV. The experi- 



mental PE spectrum (arb. units) is taken from Ref. 



30]. 
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FIG. 10: Temperature dependence of the electrical resistivity of a 239 PuCoGas single crystal after 
thermal annealing of the self-radiation damage. 
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FIG. 11: (a) Temperature dependence of the conductance curves of a contact with Rjy = 6.9Q on a 
242 PuCoGa5 crystal with reduced T c =14.5 K. The curves are vertically offset for clarity, (b) Two 
examples of low-temperature normalized conductance curves in two contacts on the same crystals 
(symbols). Lines represent the best-fit of the spectra within the d-wave 2D-BTK model. The normal- 
state resistance is Rjy = 6.9 f2 (top) and Rn = 18.8 (bottom), (c) Temperature dependence of the 
gap extracted from the fit of the conductance curves of the contact with Rn = 18.8 f2 (symbols). The 
solid line is the corresponding A(T) function calculated within the Eliashberg theory by assuming 
Qq = 6.5 meV and a coupling constant A = 2.37. 
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